
chooseBioCmirror()
5
chooseCRANmirror()
20

# if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
# 
depens<-c('GSEABase', 'GSVA', 'cancerclass', 'mixOmics', 'sparrow', 'sva' , 'ComplexHeatmap' )
for(i in 1:length(depens)){
  depen<-depens[i]
  if (!requireNamespace(depen, quietly = TRUE))  BiocManager::install(depen,update = FALSE)
}
# 
# 
# if (!requireNamespace("CoxBoost", quietly = TRUE))
#   devtools::install_github("binderh/CoxBoost")
# 
# if (!requireNamespace("fastAdaboost", quietly = TRUE))
#   devtools::install_github("souravc83/fastAdaboost")
# 
# if (!requireNamespace("Mime", quietly = TRUE))
#   devtools::install_github("l-magnificence/Mime")
# BiocManager::install(c("preprocessCore", "DESeq2"))
# 
# if (!requireNamespace("IOBR", quietly = TRUE))
#   devtools::install_github("IOBR/IOBR")
# 
# 
# devtools::install_local("Matrix")
# devtools::install_local("immunedeconv-master")
# snowfall
# BiocManager::install("snowfall")
library(snowfall)
library(devtools)
library(CoxBoost)
library(fastAdaboost)
library(Mime1)
library(IOBR)
library(dplyr)
library(immunedeconv)
library(xgboost)
# 需要准备的数据
surv = read.table(file = 'CGGA_OS.txt', sep = '\t', header = TRUE) 
colnames(surv)[1] <- "OS"
colnames(surv)[2] <- "OS.time"
surv <- surv[,-3]
rt1 = read.table("CGGA_FPKM.txt", sep = "\t")
rownames(rt1) <- rt1[,1]
colnames(rt1) <- rt1[1,]
rt1 <- rt1[-1,-1]
expr <- rt1
comgene <- intersect(colnames(expr),rownames(surv))
table(substr(comgene,14,16))
expr <- expr[,comgene]
surv <- surv[comgene,]
deg_expr <- expr %>% t() %>% as.data.frame()
surv.expr <- cbind(surv,deg_expr)
surv.expr <- cbind(row.names(surv.expr), surv.expr)
colnames(surv.expr) <- c("ID", names(surv.expr)[-1])
Dataset2 <- surv.expr
surv = read.table(file = 'TCGA-GBM.survival.tsv', sep = '\t', header = TRUE) 
surv2 = read.table(file = 'TCGA-LGG.survival.tsv', sep = '\t', header = TRUE)
merged_df <- rbind(surv, surv2)
sorted_df <- merged_df[order(merged_df[,4], decreasing = FALSE),]
surv <- sorted_df

#整理生存信息数据
surv$sample <- gsub("-",".",surv$sample)
rownames(surv) <- surv$sample
surv <- surv[,-1]
surv <- surv[,-2]
rt = read.table("GLIOMA_FPKM_MRNA2.txt", sep = "\t")
expr <- rt
colnames(expr) <- substr(colnames(expr), 1, 16)
# 替换列名中的 "-" 为 "."
colnames(expr) <- gsub("-", ".", colnames(expr))

comgene <- intersect(colnames(expr),rownames(surv))
table(substr(comgene,14,16))
expr <- expr[,comgene]
surv <- surv[comgene,]
deg_expr <- expr %>% t() %>% as.data.frame()
surv.expr <- cbind(surv,deg_expr)
surv.expr <- cbind(row.names(surv.expr), surv.expr)
colnames(surv.expr) <- c("ID", names(surv.expr)[-1])
Dataset1 <- surv.expr
# 假设Dataset1是一个数据框或数据表格
# 首先，确保Dataset1存在并且是一个数据框或数据表格
if (exists("Dataset1") && is.data.frame(Dataset1)) {
  # 互换第二列和第三列的位置
  Dataset1 <- Dataset1[, c(1, 3, 2, 4:ncol(Dataset1))]
  
  # 如果需要的话，可以将结果重新赋值给Dataset1
  # Dataset1 <- swapped_dataset
} else {
  # 如果Dataset1不存在或者不是数据框，需要处理异常情况
  print("Dataset1不存在或不是一个数据框。请检查数据。")
}
if (exists("Dataset2") && is.data.frame(Dataset2)) {
  # 互换第二列和第三列的位置
  Dataset2 <- Dataset2[, c(1, 3, 2, 4:ncol(Dataset2))]
} else {
  print("Dataset2不存在或不是一个数据框。请检查数据。")
}

list_train_vali_Data <- list(Dataset1 = Dataset1,
           Dataset2 = Dataset2)

list <- read.table("list.txt", header = F,sep = "\t", quote = "", check.names = F)
genelist <- list$V1

res <- ML.Dev.Prog.Sig(train_data = list_train_vali_Data$Dataset1,
                       list_train_vali_Data = list_train_vali_Data,
                       unicox.filter.for.candi = T,
                       unicox_p_cutoff = 0.05,
                       candidate_genes = genelist,
                       mode = 'all',
                       nodesize =5,
                       seed = 123)
custom_colors <- c("#4195C1", "#ecf0f1", "#CB5746","#3fa0c0", "#d5d9e5")
hm <- cindex_dis_all(
  res,
  color = custom_colors,       # 传递自定义颜色设置
  validate_set = "Dataset1",   # 指定测试集的名称
  order = names(list_train_vali_Data),
  width = 0.35
)
Cindex.res <- res[["Cindex.res"]]
pdf(file.path( "heatmap1.pdf"), width = cellwidth * 4 + 5, height = cellheight * nrow(Cindex.res) * 0.1)
print(hm)
dev.off()
riskscore <- res[["riskscore"]]
riskscore_1 <- riskscore[["StepCox[forward] + GBM"]]
riskscore_1_TCGA <- riskscore_1[["Dataset1"]]
riskscore_1_CGGA <- riskscore_1[["Dataset2"]]
write.table(riskscore_1_TCGA,file="TCGAriskcore.txt",sep = "\t",row.names = T,col.names = NA,quote = F)
write.table(riskscore_1_CGGA,file="CGGAriskscore.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

# cindex_dis_select(res,
#                   model="StepCox[forward] + GBM", # 指定选择模型
#                   order= names(list_train_vali_Data))
# 
# all.auc.1y <- cal_AUC_ml_res(res.by.ML.Dev.Prog.Sig = res,
#                              train_data = list_train_vali_Data[["Dataset1"]],
#                              inputmatrix.list = list_train_vali_Data,
#                              mode = 'all',
#                              AUC_time = 1,   # 一年生存
#                              auc_cal_method="KM")
# all.auc.3y <- cal_AUC_ml_res(res.by.ML.Dev.Prog.Sig = res,train_data = list_train_vali_Data[["Dataset1"]],
#                              inputmatrix.list = list_train_vali_Data,mode = 'all',AUC_time = 3,
#                              auc_cal_method="KM")
# all.auc.5y <- cal_AUC_ml_res(res.by.ML.Dev.Prog.Sig = res,train_data = list_train_vali_Data[["Dataset1"]],
#                              inputmatrix.list = list_train_vali_Data,mode = 'all',AUC_time = 5,
#                              auc_cal_method="KM")
# 
# 
# ##plot 1-year AUC predicted by all models:
# auc_dis_all(all.auc.1y,
#             dataset = names(list_train_vali_Data),
#             validate_set="Dataset1",
#             order= names(list_train_vali_Data),
#             width = 0.35,
#             year=1)
# ##Plot ROC of specific model among different datasets:
# roc_vis(all.auc.1y,
#         model_name = "StepCox[forward] + GBM",
#         dataset = names(list_train_vali_Data),
#         order= names(list_train_vali_Data),
#         anno_position=c(0.65,0.55),
#         year=1)
# #Meta-analysis of univariate cox regression for specific model
# unicox.rs.res <- cal_unicox_ml_res(res.by.ML.Dev.Prog.Sig = res,optimal.model = "StepCox[forward] + GBM",type ='categorical')
# metamodel <- cal_unicox_meta_ml_res(input = unicox.rs.res)
# meta_unicox_vis(metamodel,
#                 dataset = names(list_train_vali_Data))
# 
# ##
# rs.glioma.lgg.gbm <- cal_RS_pre.prog.sig(use_your_own_collected_sig = F,type.sig = c('LGG','GBM','Glioma'),
#                                          list_input_data = list_train_vali_Data)
# HR_com(rs.glioma.lgg.gbm,
#        res,
#        model_name="StepCox[forward] + GBM",
#        dataset=names(list_train_vali_Data),
#        type = "categorical")
# cc.glioma.lgg.gbm <- cal_cindex_pre.prog.sig(use_your_own_collected_sig = F,type.sig = c('Glioma','LGG','GBM'),
#                                              list_input_data = list_train_vali_Data)
# p2 <- cindex_comp(cc.glioma.lgg.gbm,
#             res,
#             model_name="StepCox[forward] + GBM",
#             dataset=names(list_train_vali_Data))
# pdf("多队列对比.pdf", width = 10, height = 12)
# print(p2)
# dev.off()
# 
# auc.glioma.lgg.gbm.1 <- cal_auc_pre.prog.sig(use_your_own_collected_sig = F,
#                                              type.sig = c('Glioma','LGG','GBM'),
#                                              list_input_data = list_train_vali_Data,AUC_time = 1,
#                                              auc_cal_method = 'KM')
# auc_comp(auc.glioma.lgg.gbm.1,
#          all.auc.1y,
#          model_name="StepCox[forward] + GBM",
#          dataset=names(list_train_vali_Data))
# 
# 
# ##
# for (i in c(1:2)) {
#   print(survplot[[i]]<-rs_sur(res, model_name = "StepCox[both] + plsRcox",dataset = names(list_train_vali_Data)[i],
#                               #color=c("blue","green"),
#                               median.line = "hv",
#                               cutoff = 0.5,
#                               conf.int = T,
#                               xlab="Day",pval.coord=c(1000,0.9)))}

res.feature.all <- ML.Corefeature.Prog.Screen(InputMatrix = list_train_vali_Data$Dataset1,
                                              candidate_genes = genelist,
                                              mode = "all",nodesize =5,seed = 5201314 )
write.table(res.feature.all,"所有算法筛选对象.txt", col.names = T, row.names = F, sep = "\t", quote = F)


core_feature_select(res.feature.all) ##空的
##实现上一步骤 Upset图
library(UpSetR)
library(ggplot2)

# 假设 res.feature.all 是一个包含方法和所选特征的数据框
# 例如：
# res.feature.all <- data.frame(
#   method = c("method1", "method1", "method2", "method2"),
#   selected.fea = c("gene1", "gene2", "gene3", "gene4")
# )

# 设置颜色和比例
col <- c("#E18727", "#B09C85", "#ADB6B6", "#B09C85")
mb.ratio <- c(0.6, 0.4)

# 创建核心特征列表
core_feature_list <- list()
for (i in unique(res.feature.all$method)) {
  core_feature_list[[i]] <- res.feature.all[res.feature.all$method == i, "selected.fea"]
}

# 创建 UpSet 图
p1 <- upset(
  fromList(core_feature_list), 
  sets = names(core_feature_list), 
  order.by = "freq", 
  nintersects = NA, 
  mb.ratio = mb.ratio, 
  keep.order = TRUE, 
  mainbar.y.label = "Shared gene number", 
  sets.x.label = "Total gene number", 
  point.size = 2, 
  line.size = 1, 
  sets.bar.color = col[1], 
  main.bar.color = col[2], 
  matrix.color = col[3], 
  shade.color = col[4]
)

# 显示 UpSet 图
print(p1)
pdf("core_feature_upset_plot.pdf", width = 10, height = 6)
print(p1)
dev.off()

core_feature_rank(res.feature.all, top=20)
pdf(file.path( "heatmap1.pdf"), width = cellwidth * 4 + 5, height = cellheight * nrow(Cindex.res) * 0.1)

test.sig <- signature_score(res, list_train_vali_Data)
save.image(file = "current_workspace.Rdata")